##########################################################################################

library(data.table)
library(optparse)
library(parallel)
library(dplyr)
library(RColorBrewer)
library(ggplot2)
library(ggrepel)

##########################################################################################

option_list <- list(
    make_option(c("--input_file"), type = "character"),
    make_option(c("--gtf_file"), type = "character"),
    make_option(c("--q_t"), type = "character"),
    make_option(c("--foldchange_t"), type = "character"),
    make_option(c("--out_path"), type = "character")
)

if(1!=1){
    
    input_file <- "~/20220915_gastric_multiple/rna_combine/analysis/images/DiffGene/DiffGene.tsv"
    out_path <- "~/20220915_gastric_multiple/rna_combine/analysis/images/DiffGene"
    q_t <- 0.05
    foldchange_t <- 2
    gtf_file <- "~/ref/GTF/gencode.v19.ensg_genename.txt"

}

parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

input_file <- opt$input_file
out_path <- opt$out_path
gtf_file <- opt$gtf_file
q_t <- as.numeric(opt$q_t)
foldchange_t <- as.numeric(opt$foldchange_t)

dir.create(out_path , recursive = T)

##########################################################################################

dat_diff <- data.frame(fread(input_file))
dat_gtf <- data.frame(fread(gtf_file , header = F))
colnames(dat_gtf) <- c("gene_id" , "Hugo_Symbol")
dat_diff <- merge( dat_diff , dat_gtf , by = "gene_id" )

##########################################################################################

class_type <- c( "Normal" , "SG" , "AG" , "IM" , "DYS" , "IGC" , "DGC")

image_name <- paste0( out_path , "/DiffGene.HugoSymbol.tsv" )
write.table( dat_diff , image_name , row.names = F , quote = F , sep = "\t" )

##########################################################################################
diffplot <- function(use_dat = use_dat , image_name = image_name){

    #对原数据进行处理
    use_dat$log10_q_Value <- -log10(use_dat$padj)
    use_dat$log_FC <- use_dat$log2FoldChange
    use_dat$gene <- use_dat$Hugo_Symbol

    data <- use_dat[,c('gene','log_FC','log10_q_Value')]
    colnames(data) <- c('gene','log_FC','-log10_P_Value')

    #设置阈值
    logFC_cutoff <- log2(foldchange_t)
    log10_P_Value_cutoff <- -log10(q_t)

    options(ggrepel.max.overlaps = 20)

    plot1 <- ggplot(data = data,aes(x = log_FC,y = `-log10_P_Value`))+
      geom_point(data = subset(data,abs(log_FC)<logFC_cutoff),
                 aes(size = abs(log_FC)),col = 'gray',alpha = 0.4)+
      geom_point(data = subset(data,abs(`-log10_P_Value`)<log10_P_Value_cutoff & abs(log_FC)>logFC_cutoff),
                 aes(size = abs(log_FC)),col = 'gray',alpha = 0.4)+
      geom_point(data = subset(data,abs(`-log10_P_Value`)>log10_P_Value_cutoff & log_FC>logFC_cutoff),
                 aes(size = abs(log_FC)),col = 'red',alpha = 0.4)+
      geom_point(data = subset(data,abs(`-log10_P_Value`)>log10_P_Value_cutoff & log_FC< -logFC_cutoff),
                 aes(size = abs(log_FC)),col = 'darkgreen',alpha = 0.4)+
      theme_bw()+
      theme(legend.title = element_blank(),
            panel.grid.major = element_blank(),
            panel.grid.minor = element_blank(),
            legend.position = 'none',
            axis.line = element_line(colour = "black"))+
      labs(x='log2(fold change)',y='-log10(adjusted p-value)')+
      geom_vline(xintercept = c(-logFC_cutoff,logFC_cutoff),lty = 3,col = 'black',lwd = 0.4)+
      geom_hline(yintercept = log10_P_Value_cutoff,lty = 3,col = 'black',lwd = 0.4) +
      geom_text_repel(data = subset(data,abs(`-log10_P_Value`)>log10_P_Value_cutoff & abs(log_FC)>logFC_cutoff),
                            aes(label = gene),size = 3,col = 'black')
    ggsave( image_name , plot1 , width = 10 )
}

for( i in 1:(length(class_type)-1) ){
    for( j in (i+1):length(class_type) ){
        class_1 <- class_type[i]
        class_2 <- class_type[j]
        use_dat <- subset( dat_diff , class2==class_2 & class1==class_1)
        image_name <- paste0( out_path , "/" , class_1 , "_" , class_2 , ".DiffGene.pdf" )
        diffplot(use_dat = use_dat , image_name = image_name)
    }   
}
